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Abstract 



We use lattice Boltzmann simulations to investigate the formation of arrested struc- 
tures upon demixing of a binary solvent containing neutrally wetting colloidal particles. 
Previous simulations for symmetric fluid quenches pointed to the formation of 'bijels': 
bicontinuous interfacially jammed emulsion gels. These should be created when a glassy 
monolayer of particles forms at the fluid-fluid interface, arresting further demixing, and 
rigidifying the structure. Experimental work has broadly confirmed this scenario, but 
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shows that bijels can also be formed in volumetrically asymmetric quenches. Here we 
present new simulation results for such quenches, compare these to the symmetric case, 
and find a crossover to an arrested droplet phase at strong asymmetry. We then make 
extensive new analyses of the post-arrest dynamics in our simulated bijel and droplet 
S^ structures, on time scales comparable to the Brownian time for colloid motion. Our 

results suggest that, on these intermediate time scales, the effective activation barrier 
to ejection of particles from the fluid- fluid interface is smaller by at least two orders of 
magnitude than the corresponding barrier for an isolated particle on a flat interface. 

1 Introduction 

The tendency of colloidal particles to sequester at fluid-fluid interfaces has been known for 
over a century, originating with the work of Pickering [lj and Ramsden [2] who used colloids 



to stabilize droplet emulsions. Recently, understanding of the formation of such particle- 
stabilized emulsions has advanced considerably [3]. These advances, combined with the in- 
creasing availability of colloidal nanoparticles (in some cases with tunable surface chemistry) 
[3J, H] , has led to strong interest in using liquid-liquid interfaces to direct the self-assembly 
of such nanoparticles [U G3 [6]. In most such studies the liquid-liquid interface is created 
by agitation or sonication of the demixed fluids [3J, H] . This generally leads to formation of 
simply connected (droplet) structures, although these can exhibit frozen aspherical shapes 
[3I| (a phenomenon also seen in particle-stabilized gas bubbles 0). In such cases, the 
particle layer has solidified, imparting mechanical rigidity to each droplet. 

Recent simulation work [10] addressed a new route to directed self assembly of nanocol- 
loids, creating fluid-bicontinuous structures. If the interface solidifies in such a case it should 
impart macroscopic rigidity in three dimensions [TOl ITT]. The resulting 'bijel' (bicontinuous 
interfacially jammed emulsion gel), comprising a solid matrix permeated by a pair of bi- 
continuous fluids, could have potential applications as a 'membrane contactor' for catalytic 
applications pjjl [121 [13] . I n the published simulation protocol [10], one starts with a sample 
in the single phase region of a symmetric (50:50 by volume) binary fluid, in which are sus- 
pended colloidal particles having equal affinity for the two fluids. There follows a quench in 
which the fluids demix, sweeping up particles to the interface; since the solid-liquid interfacial 
tensions with the two solvents are equal, these adopt a 90° contact angle (neutral wetting). 
As the interface reduces its area by the usual coarsening process, the colloids become jammed 
into close proximity, and coarsening is dramatically curtailed [10]. An alternative mesoscopic 
simulation method, dissipative particle dynamics, has very recently been used [14] to confirm 
several findings of the earlier lattice Boltzmann studies [10] ■ 

Although these simulations cannot be run for long enough to determine the ultimate fate 
of the resulting structure, arguments were given [10] for its permanent arrest. Specifically, 
there is an energy barrier cue, where e = ana 2 , to detachment of a particle of radius a from 
a fluid-fluid interface of tension o. We may write e/fc^T = (a/ao) 2 where a 2 , = hsT/ira; 
then for T = 300 K and typical a of order 0.01 Nm" 1 or larger, a is 0.4 nm or less. The 
geometry-dependent parameter a is of order unity for a single particle on a flat surface, 
and arguably should be similar in a dense particle layer; if so, ae/ksT > 10 even for a 
particle of 1 nm radius, and thermally activated detachment can be safely neglected for, say, 
a > 3 nm. On the other hand, due to the complicated energy landscape involved, we cannot 
quantitatively estimate a for a crowded monolayer layer with any precision. The possibility 
remains that such a layer might (unlike a single adsorbed particle) have pathways allowing 
sequential particle expulsions while continuously decreasing the fluid-fluid interfacial area. 
This would correspond to a — 0. 

Encouragingly, fully arrested bijels have recently been created experimentally [15J. As 
predicted [JU] , they are soft solids and exhibit a yield stress. The first examples were prepared 
in thin slabs and had rather ill-formed fluffy colloidal multilayers at the fluid-fluid interface 
[8]. (Similar slab-based work on polymer blends has also been reported [16J.) However, 
improved quenching technique has recently yielded bulk three-dimensional samples with a 
clean colloidal monolayer instead [15]. These high-quality bijels nonetheless deviate from the 



idealized protocol simulated by Stratford et al. [10] in several respects. First, the colloidal 
particles used - although repulsive at moderate separations - might on close approach fall 
into a primary flocculation minimum creating quasi-permanent bonds. They are also much 
larger (a = 290nm) than the simulated nanocolloids [10], reducing the possible role of Brow- 
nian motion. The particles do not have exactly neutral wetting, and thus do not sit quite 
symmetrically across the fluid-fluid interface. Finally, the binary fluid mixture is itself not 
symmetric: the fluids have different viscosities, a somewhat asymmetric phase diagram, and 
are quenched to form unequal volume fractions of the two phases [Bl [T2] . 

In this paper we present in Section [3] new simulation studies of the crossover from bicon- 
tinuous to droplet morphologies, focussing for simplicity on the case of quench asymmetry. 
This is preceded in Section [2] by a brief account of the simulation methodology. We then 
quantify in Section H] the residual dynamics at intermediate times (of order the Brownian 
time for a particle) of both the fluid domains and the colloidal particles in our bijel and 
arrested droplet morphologies. This includes a study of the temperature dependence of par- 
ticle ejection rates, from which we deduce that, on these time scales at least, the barrier 
parameter a is surprisingly small. Our conclusions are in Section [51 

2 Simulation Methods 

We use the lattice Boltzmann method for a binary fluid incorporating spherical solid particles 
|17j . modifying a standard 'bounce-back on links' method [18], to allow for the presence of 
a binary solvent p21 [19]. The two solid-fluid interfacial tensions are exactly equal, with 
the interfacial thermodynamics implemented as reported previously [T7_]. We choose for 
simplicity a pair of fluids with equal density p and viscosity rj. Their phase diagram is also 
symmetric, and described by the free energy functional [20] 



dv(A^j 2 /2 + B^/A + K(Vij) 2 /2) (V 



where the order parameter if) describes the fluid composition, and the model parameters 
A < 0, B, and k control the fluid-fluid interfacial tension a and thickness £ [20] . The if) = 
isosurface defines the position of the fluid-fluid interface. 

The binary fluid is initialised to be well mixed and at rest, with a small amplitude 
random noise added to the if) field to induce spinodal decomposition. The colloids are 
initially positioned at rest randomly throughout the system. We choose a deep quench in 
the sense that thermal fluctuations of if) remain negligible: without particles, the phase 
separation proceeds purely by amplification of the initial noise. Time-dependent, thermal 
noise is however fully included in the description of fluid momentum, using a method reported 
previously [21]. As a result, the colloids undergo realistic, many-body Brownian motion. 

Unless otherwise stated below, our parameters are set as follows. First, we choose —A = 
B = 0.002, and k = 0.0014, giving an interfacial thickness of £ = 1.14 and tension a = 
1.58 x 10 -3 ; we set the (scalable) fluid density p — 1, and choose viscosity rj = 0.1. Here and 
below, all physical variables are expressed in lattice units [20] (except where SI units are given 



explicitly). For numerical efficiency, the colloid radius a is made as small as is compatible 
with acceptably accurate simulation [TQI [22]; we choose a = 2.3, so that e = 0.026. Our 
noise setting is k B T = 2.13 x 10~ 5 , so that e/k B T = 1230. We can choose to interpret 
these parameters as representing a short-chain hydrocarbon/water or hydrocarbon/alcohol 
mixture, with laboratory parameters p = 1000 kg m -3 , t] = 9.3 x 10 -4 N m _2 s (close to 
that of water at ~ 300K) and o = 6.1 x 10~ 2 N m _1 . This corresponds to a physical particle 
radius of a = 5.1 nm. 

In the absence of colloidal particles, the only characteristic length and time scales as- 
sociated with the physics of hydrodynamic coarsening for binary fluids are [20j [23] L = 
rj 2 /{per) and t = rj 3 /(pa 2 ), which for the fluid parameters just chosen are L « 14 nm 
and t — 0.22 ns. Computing the same quantities in lattice units (L = 6.33, t = 401) 
allows length and time scales to be matched to experiment, in principle. In practice we fully 
match e/ksT = 1230 and a/L = 0.363; together these also ensure matching of r^/to, where 
tb = 67Tr]a^/kBT is the Brownian time for a free colloid to diffuse its own radius. How- 
ever, not all dimensionless control parameters of potential relevance to the problem can be 
fully matched. For instance, the particle-scale Reynolds number, Re = (dL/dt)pa/r], which 
characterises the relative importance of fluid inertia to viscosity, cannot be made as small 
as the true physical value, but can be made small compared to unity, which we take to be 
sufficient [22]- Here L(t) is the domain-scale correlation length of the demixing fluid, which 
we conventionally define as: 

w-Hm& (2) 

where S(k, t) = (\ip(k, t)\ 2 ) is the equal time structure factor and (.) denotes an average over 
a shell in k-space at fixed k = |k|. (For definiteness, we set ip = in the interior of the 
particles.) 

As viewed by the lattice fluid, our colloids are spherical only on time average. To take 
account of this, a calibration is performed to find their hydrodynamic radius a^ [18]. This 
is the radius of the sphere which exhibits the same mean Stokes friction coefficient Qirrjah 
as the lattice colloid. We have chosen a combination of viscosity and particle size such that 
the actual radius (defined by the stencil that creates the bounceback links) and a^ coincide: 
a = ah = 2.3. Another effect of the discretization is that hydrodynamic interactions between 
two colloidal particles are under-represented when their surface-to-surface separation is small 
on the lattice scale. This allows unphysical overlap of particles, risking numerical instability. 
One can rectify this by adding lubrication forces by hand [T7J [18], at the cost of much 
increased run-times whenever there are multiple lubrication contacts between particles. To 
avoid that problem, a soft-core thermodynamic repulsion can be added so that particles 
maintain a minimum separation 2a? somewhat larger than 2a,h [TU]. (Such a repulsion is 
often present physically, e.g., for charge-stabilized colloids.) 

For the current work, a number of exploratory simulations were performed on a D3Q15 
lattice of volume A 3 = 64 3 with periodic boundary conditions, using exactly the same com- 
bination of lubrication and soft-core forces as was described by Stratford et al. [10]. (DNQM 
signfies an N-dimensional lattice with M discrete velocities.) To obtain the production- run 



datasets reported below, we used a larger D3Q19 lattice (A 3 = 128 3 ), and deployed an im- 
proved soft-core approximation to the hard-core repulsion [24]. Despite a somewhat smaller 
ax, this has better stability properties. Also, the explicit lubrication forces, which in the 
exploratory work were found to barely influence the arrested-state properties, were switched 
off in these production runs. In what follows, the colloid volume fraction <p is always defined 
in terms of a^: = Amra^/3 where n is the colloid number density. 

3 Coarsening and Arrest Dynamics 

If the binary fluid mixture in which particles reside is made sufficiently asymmetric, one 
expects a crossover from bicontinuous to droplet-like structures. This asymmetry could be 
brought about via viscosity (forming droplets of the less viscous fluid [23]); via wetting 
(forming droplets of the less wetting fluid); via phase-diagram asymmetry (forming droplets 
of the purer of the two phases); or via quench asymmetry (forming droplets of the minority 
phase). In experimental practice, these tendencies might be played off against one another 
[8] and each mechanism could have its own subtle morphological characteristics. 

In the present work, we focus purely on quench asymmetry, and vary the mean initial 
order parameter ip . With the symmetric free energy ([T]), values of ip = 0,0.1,0.2,0.3,0.4 
give respectively phase volume ratios 50:50, 55:45, 60:40, 65:35 and 70:30. (Because of 
the remaining symmetries, the sign of ip can be reversed without physical consequence.) 
In all cases we perform a deep quench so that the mode of phase separation is spinodal 
decomposition, as opposed to nucleation and growth |26j. 

3.1 Domain Morphology 

In the spinodal demixing of a colloid-free binary fluid, there is a threshold ipp for the loss 
of bicontinuity which has been estimated theoretically, experimentally, and by simulation as 
ip p — 0.44±0.04 (phase volume 0.28±0.02 |26j). For ip moderately beyond ip p , on quenching 
into the spinodal regime, an initially bicontinuous domain pattern is formed by diffusion, 
but once the late (hydrodynamic) stages are entered, this breaks into droplets. 

When colloids are added, we do not expect these to perturb strongly the coarsening 
dynamics until their density on the fluid-fluid interface is high enough to cause jamming. 
Our simulations place this time well into the hydrodynamic phase of the coarsening process, 
so we expect a similar ip p in the presence or absence of particles. This is indeed observed. 
In Figfj] we present snapshot configurations at t = 5 x 10 5 lattice units (corresponding to 
t = 275 ns), for ip = 0.1,0.2,0.3, with particle volume fraction <fi = 0.2. In all these 
configurations the structure has substantially stopped evolving, although there is a residual 
slow dynamics, discussed separately below in Section HI All these cases remain bicontinuous 
throughout the simulation. In contrast, Fig|2J shows a similar quench but with ip = 0.4. By 
the time of arrest, the system has now depercolated into an array of deformed droplets, of 
essentially frozen shape. The largest droplets present are smaller than the system size (note 



that only one eighth of the system is shown). Thus in our simulations with colloids present, 
0.3 < ip p < 0.4, suggesting a modest decrease, if any, in ip p from the colloid-free case. 




Figure 1: Snapshot late time configurations for = 0.2 and ip = 0.1, 0.2, 0.3 (top to bottom); 
lattice size A = 128, cropped to A = 64. The right column shows the fluid-fluid interface 
plus particles rendered as spheres of the appropriate hydrodynamic radius. The left column 
shows the same image with particles made transparent to aid visualization of the domain 
topology. Parameter settings as given in methods section. 



3.2 Domain Growth Kinetics 

The time evolution of L(t) as defined by Eqj2]is shown for a symmetric quench with colloid 
volume fraction <fi = 0.2 in Figj3l The evolution is similar to that reported previously [10] 
for the same <ft an d lattice size. (The L values are somewhat larger, as expected for the 
slightly reduced thermodynamic radius, a?, in the present work.) The run time shown 
(10 6 ) extends 50% beyond any previously published; however with our chosen parameter 
mapping it still represents less than 1 fis in laboratory time. We thus attain, but do not 
significantly exceed, the colloidal Brownian time tb — 10 6 . As is clear from the plots, 
during the newly investigated time window, L(t) continues to increase at an ever slowing 
rate. Eventual saturation is not seen, but nor can it be ruled out. In Figj3] we also show 




Figure 2: Snapshot late time configuration for <fi = 0.2 and ip = 0.4; lattice size A = 128, 
cropped to A = 64; interface and particles represented as in Fig. [Q 
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Figure 3: Time evolution of domain size L(t) for various quenches; lattice size A = 128. 
Symmetric quench ip = (topmost curve) with colloid volume fraction = 0.2; asymmetric 
quenches of ip = 0.4 with <fi = 0.2 (middle curve) and 4> — 0.25 (lowest curve). Parameter 
settings as given in methods section. 



data for two asymmetric quenches with ip = 0.4 and 4> — 0.2, 0.25. Both lead to droplet 
morphologies resembling that in FigO The slow residual dynamics seen in the symmetric 
quench is still present, at least for <fi = 0.2, despite the loss of fluid bicontinuity. The run 
with <p = 0.2, ip = 0.4 was further extended to t = 1.4 x 10 6 , resulting in a continued gradual 
increase of domain size up to L = 37 (data not shown). Note that saturating trend in the 
datasets of Figj3] is not well represented by any power law; replotting the data on a log-log 
representation (not shown) still yields significant curvature throughout the range shown. 

4 Dynamics at Intermediate Times 

In this section we present various data analyses intended to shed more light on the residual 
dynamics responsible for the continuing increase of L(t) for 2 x 10 5 < t < 10 6 as seen in 
FigEJ Since we do not achieve t ^> tb, we refer to this as the 'intermediate' time regime, 
leaving open the question of what happens ultimately. 

We can think of three candidate explanations for the residual dynamics, as follows: (i) 
Formation of semi-crystalline rafts at the interface, allowing closer packing of particles; 
(ii) Continuous ejection of particles from the interface due to a very small, possibly even 
vanishing, value of the geometric barrier parameter a (see introduction) for a crowded layer; 
(iii) Numerical artefact, e.g. due to inadequate resolution of fluid interfaces, causing the 
effective a in the simulations to be less than it should be. 

Of these, the first process should be ultimately self-limiting; and it can be suppressed 
by choosing bidisperse particles. Doing so reduces dL/dt at late times, but not to zero [10J; 
instead of slow crystallization, there is a gradual segregation by size which might prove even 
slower. We therefore do not pursue the bidisperse case in the work presented here. The 
second explanation is interesting, but somewhat at odds with the experimental observation 
that bijels can remain stable for weeks [8j[l5]. Recall however that the experimental systems 
involve much larger particles [TH] which might also be stabilized, once jammed on the inter- 
face, by attractive interparticle forces. Moreover, as we discuss below in Section I4.5[ a could 
effectively be time dependent, achieving large values only for t^> tb- The final explanation, 
numerical artefact, cannot be completely excluded. In our simulations, a is not much larger 
than £, the interfacial width; failure to separate these lengths adequately might exaggerate 
the effect of diffusive particle displacements in overcoming the barrier to detachment. When 
structural motifs (ripples and cylinders coated with colloids) were simulated at higher res- 
olution than attempted here |10j . they did not show much particle ejection. This result is 
however not conclusive for bijels since a might depend on interfacial shape; moreover, these 
runs could only probe t <C tb- To settle the issue definitely would require simulations of a 
bulk bijel with a ~^> 2.3. This remains out of reach with current resources; note that a factor 
two radius increase requires nearly a tenfold increase in computational effort. The results of 
FigHUdo appear to settle one point however. Since the slow increase in L(t) continues even 
in the disconnected droplet phase, failure to adequately resolve thin fluid necks at domain 
pinchoff [10J is not the main cause of the residual coarsening. 



4.1 Interfacial Ordering 

To quantify the ordering effect, we show in FigHJthe intermediate-time radial distribution 
function g(r) for the colloids, comparing the case of symmetric (bicontinuous) and asym- 
metric (droplet structure, ip = 0.4) quenches. Note that g{r) is defined as usual for the 3D 
bulk sample: we make no attempt to restrict examination to the interfacial environment. 
However, sharp features in g{r) at length scales up to a few a T are presumably a probe of 
interfacial ordering since only on the interface are particles in close proximity. During the 
period where L is evolving slowly, there is significant sharpening of the peaks in g(r) without 
significant changes in their position. This is consistent with the slow formation of crystalline 
rafts, and can be quantified by the time evolution of the height of the first peak (FigfJJ). Also 
shown is g{r) in a symmetric quench after thermal noise was switched off mid-run. The effect 
of this is to allow a period of downhill relaxation to a local minimum of the (interfacial plus 
interparticle) energy. This further sharpens the peaks with very little change in position. 
This suggests that the particle layer was effectively trapped in a jammed state, very close to 
a local metastable minimum, even prior to switching off the noise. 

4.2 Particle Ejection 

To quantify the ejection of particles we plot in Figj5] the time evolution of the number of 
'free' particles, Nf(t). A free particle is defined as one which is not in contact with the fluid- 
fluid interface. Interfacial contact is identified by smallness of the local order parameter 
ip at nodes contacting the particle surface. (However, in a crowded layer, many of these 
nodes are occupied by other particles. Because of this and other discretization effects, it 
is not practical within studies at this resolution to classify the non-free particles by the 
strength of their binding to the interface.) These results, for both symmetric and droplet 
morphologies, show that local crystallization alone is not enough to explain the residual 
dynamics at intermediate times. 

For the symmetric quench, the number of free particles increases monotonically through- 
out the slow dynamics regime. (There is statistical noise associated with particles that 
transiently satisfy the criterion for inclusion in Nf without really escaping the interface.) 
The same is seen for the droplet quench, where we also find that Nf is notably higher. In 
this case, the free particles are more numerous in the continuous phase than in the droplets, 
suggesting preferential detachment from the exterior surface of droplets. This is physically 
reasonable for a crowded curved layer in which the interparticle forces will tend to displace 
colloid centres towards the exterior. 

In Figj5]we also show the effect on Nf(t) of switching off thermal noise at t — 4 x 10 5 in 
a symmetric quench. For a significant period after switchoff, there is no continued increase 
of Nf(t). However there is some indication of ejection resuming at very late times. The 
L(t) curve for the same system is also shown; notably, although switching off the noise does 
reduce the rate of domain growth, this reduction is very modest and the coarsening rate 
remains finite throughout the period of non-ejection. 
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Figure 4: Late time g(r) for colloidal particles; the data shown is time averaged over the 
interval 9.5 x 10 5 < t < 10 6 . (a) ip = 0.4; (b) ip = 0; (c) ip = 0, with thermal noise 
switched off at t — 4 x 10 5 ; (d) time evolution of the height of the first peak for conditions 
(a)-(c) (bottom to top), each data point averaged over a 5 x 10 4 timestep bin. Parameter 
settings as given in methods section. 
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Figure 5: (a) Time evolution of the number of free particles Nf(t) for droplets (ip = 0.4 
upper curve) and a symmetric quench (lower curve), each with A = 128,0 = 0.2. (b) The 
same data for the symmetric quench (upper curve) compared with the case where noise is 
switched off after t = 4 x 10 5 (lower curve), (c): L(t) data for the two systems in (b), 
upper curve with noise, lower with noise switchoff at t — 4 x 10 5 , inset zoom in on late 
time regime. Parameter settings as given in methods section. Note that our definition of Nf 
entails Nf = in a system prior to demixing (in which ip ~ everywhere). However, for 
all the datasets shown sharp interfaces are present once t > 2 x 10 4 , so that Nf(t) has the 
stated meaning for later times than this. 



11 



4.3 Particle Dynamics 

We show in Figj6] the RMSD (root-mean-square displacement) of the colloidal particles, 
R(t, t w ), between 'waiting times' t w and later times t w + 1. For both symmetric and droplet 
morphologies, the RMSD at any fixed t decreases monotonically with t w , as expected in a 
system that is continuously slowing down. There is evidence of diffusive behavior at large 
t w , although R values remain small on the scale of the particle radius a. Additionally 
there are signs of upward curvature at large t at the higher values of t w . We also show in 
Figj6] a symmetric quench where thermal noise is switched off mid-run. Upon switch-off, 
the behavior at modest t changes from diffusive (R ~ lO -3 ^ 1 / 2 ) to ballistic, but with a 
much smaller magnitude: R ~ 10 — 5 t x . The coefficient here is comparable to dL/dt in the 
intermediate time regime, suggesting that the residual motion may stem from advection of 
particles by the still coarsening fluid-fluid interface. At larger times, t > 10 4 , this ballistic 
behavior crosses over to a lower slope on the log plot, with R remaining slightly below the 
values found when the noise was left in place throughout. 

4.4 Droplet Dynamics 

We show in FigJ7|the time evolution of a single droplet as it occurs within our simulation for 
ip = 0.4. The earliest time shown is soon after the droplet adopts its final topology. It then 
moves towards a spherical structure which is however not reached, instead forming a faceted 
drop reminiscent of the colloid-armored gas bubbles reported elsewhere [9]. (Arrested non- 
spherical emulsion droplets have also been seen [TJ [8] , where however the interface remained 
amorphous.) On the facets, a clear tendency towards local crystallinity is seen. For this 
particular droplet (~ 120 particles) there is no particle ejection after formation, and little 
discernible change of any kind after t~5x 10 5 . 

4.5 Residual Dynamics: Discussion of Mechanism 

The data of Sections 14. 1J4.2I show that the slow dynamics do not solely involve consolidation 
of particles into crystalline regions, but also involve particle ejection. On the other hand, 
suppressing ejection (at least temporarily) by switching off Brownian motion entirely does 
not prevent the residual dynamics either. There seems to be some interplay between different 
mechanisms at work here. The RMSD data of Section I4T51 seem broadly consistent with the 
idea that both thermal and athermal mechanisms contribute. 

Recall that, despite the very long run times involved, the time window under study barely 
achieves the Brownian time scale tb — 10 6 . The residual dynamics could result from the 
fact that the particles are gathered onto the interface, which then arrests, on a timescale 
short compared to t#. Such particles might not achieve local thermal equilibrium, with 
respect to the interparticle potential and interfacial forces, until well after the initial arrest 
event. Continued ejection of particles then remains possible up to t ~ tb, even if a is 
not small, so long as these ejections involve a minority of particles that have only tenuous 
contact with the fluid interface itself, and thus have relatively low barriers (a few k B T) to 
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Figure 6: Log-log plots of colloid RMSD: (a,b) for t^/10 4 = 2 n , n = 0,1, 2, ...6 (upper to 
lower curves) with (a) ip = 0.4; (b) ip = 0. In (c) ip = 0, curves for t^/10 5 = 4, 6, 8 
are compared with and without thermal noise switched off at t — 4 x 10 5 . Straight lines in 
(a,b) have slope 1/2 representing diffusive behavior. Parameter settings as given in methods 
section. 
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Figure 7: Time evolution of droplet shape. Snapshots at t/10 5 values as follows (a):l, (b):2, 
(c):3, (d):4, (e):5, (f ):10, (g):14. The number of attached particles does not change during 
this sequence. 



removal. Thereafter, ejections should become ever rarer and eventually cease. In effect, 
this scenario amounts to an a value which increases over time as the population of weakly 
bonded particles is depleted. A broad distribution of barrier heights could give rise to very 
slow 'aging' dynamics on timescales tb and beyond [27]. Since t b itself rapidly becomes large 
as the colloid radius is increased, such dynamics might well be detectable on experimental 
time scales. 



4.6 Direct Estimation of Activation Barrier 

In principle, any activation barrier Ea to particle ejection can be measured simply by ob- 
serving the ejection rate r (per adsorbed interfacial particle) as a function of temperature 
and fitting to the Arrhenius form, ln(r/ro) = —EaJ^bE + constant. In practice however, 
one expects this expression for the rate to become reliable only when t ^> Tq 1 , so that the 
escape process is sampled on time scales large compared to the "attempt frequency" r for 
the barrier crossing process. The temperature dependence of r is discussed below, but the 
natural first assumption is to set tq — t^ 1 oc ksT. Note that a reduction in diffusivity 
arising from crowding effects could decrease the attempt frequency and require even longer 
time scales to be probed; but this rate reduction would be temperature- independent, and not 
directly visible in Ea- Also invisible in Ea would be any entropic, as opposed to energetic, 
barrier to particle ejection. 

In an attempt to directly estimate Ea, we have repeated our simulations for the symmetric 
quench (0 = 0.2) for seven higher temperatures, creating a temperature span of one decade. 
The Brownian time at the highest of these temperatures is tb — 10 5 , allowing us to probe 



14 




1.2e+06 



Figure 8: Log linear plot of surviving number Nt of trapped particles against time for tem- 
peratures (from top to bottom at extreme right) 10 5 k B T = 2.13, 3f, 4.5, 6, 8, llf, 12.5, 15, 20. 
Those marked f are averages over two runs; the variation in late-time slope for repeat runs 
with the same parameters is about ±10%. 
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Figure 9: Plot of lnr vs 1/kT for the four highest temperatures. The linear least-squares fit 
is shown. 
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about one decade beyond t b . In this simulation particle ejection is rather frequent, so that 
after 10r^ more than 20 percent of the particles initially adsorbed to the interface have been 
ejected. However, this temperature (2 x 10~ 4 ) is at the extreme upper end of the acceptably 
accurate range for LB [21] . Backing off to more conservative values creates a more reliable 
simulation, but with poorer statistics for the ejection process and (crucially) a longer tb- 

In Figj8]we show for all temperatures the quantity lniVV(t) plotted against time. Here 
Nt is the number of particles trapped on the interface; this is defined by N? + Nf = N, 
with N the total particle number in the simulation. For a simple activated process we 
expect a series of straight-line plots with slope —r(T). (For a local process, there should be 
no dependence of r on L, whose late-time variation is anyway modest for all these runs.) 
Roughly constant slopes are indeed apparent at late times for the higher temperatures, but in 
all cases significant curvature is present for t/r B < 3.5. To obtain reliable rate estimates we 
therefore excluded data in that early-time window. Only four substantive datasets remain; 
these were fitted by least-squares on the log-linear representation of FigJS] to find the rates 
r(T) for the four highest temperatures. 

The standard Arrhenius analysis would assume that the pre-exponential rate factor r$ 
is diffusive and hence linear in temperature. With this assumption one can estimate Ea 
by fitting the four r(T) values to a straight line on a plot of \n(r/k B T) vs l/k B T. The fit 
however is quite poor; the resulting estimate is Ea = 7 x 10~ 5 ± 2 x 10~ 5 . Interestingly, a 
much better fit is obtained by assuming r(T) is temperature- independent, as might happen 
if the barrier-crossing attempt rate were fixed by the intrinsic coarsening dynamics of the 
interface (which have velocity scale cr/rj) rather than diffusion. This fit is shown in Figj9j 
and gives E A = 2.2 x 1(T 4 ± 1.6 x 1(T 5 (with r = 5.5 x 1(T 7 ± 7 x 1(T 8 ). 

Inclusion of results from lower temperatures (necessitating use of data from t/r B < 3.5) 
would lead to systematically smaller estimates for Ea- Indeed, using the data for t > 6 x 10 5 
(say) from all datasets gives an Ea value whose difference from zero is not statistically 
significant. We conclude that any activation energy Ea measurable on the Brownian time 
scale is not much more than about 2 x 10 -4 ; and possibly much less. Thus for the simulations 
reported in previous sections, which have k B T = 2.13 x 10~ 5 , we find Ea < 10k B T. Equating 
this to at and noting that e/k B T = 1230, gives a < 8 x 10~ 3 . For a dimensionless quantity 
that was expected to be of order unity, this is suspiciously close to zero; clearly a = cannot 
be ruled out. 

As suggested above, the apparent smallness of a might stem from a penumbra of weakly- 
attached particles, so that the ejection rate in this 'late-intermediate' time window (3.5t b < 
t < IOtb) remains controlled by atypically low barriers. This would entail eventual upward 
curvature to the plots in FigJHl There is no sign of this, but the vertical scale is quite 
expanded (i.e. no more than 15 percent of particles are ejected during any of our fitted 
straight-line windows), so neither is it excluded. 



16 



5 Conclusions 

In this work we studied the crossover from bijel to frozen droplet structures on varying the 
volumes of the two demixing phases in a binary fluid system containing neutrally wetting 
nanocolloids. The depercolation threshold for ip p was at the low end of expectations based 
on previous work on binary fluids without particles. We then compared the domain growth 
kinetics for asymmetric and symmetric quenches; both show residual coarsening dynamics 
at times following the jamming transition for interfacial particles. This dynamics appears to 
involve an interplay of particle ejection and particle rearrangement to form locally crystalline 
packings; neither appears solely responsible for the residual interfacial evolution. This con- 
clusion was reached by studying the time evolution of the colloid pair distribution function, 
the number of unattached particles, and the root mean squared displacement of the col- 
loids. Studying the time evolution of a single droplet within the quenched structure shows 
relaxation towards an aspherical, faceted shape showing significant local crystallinity 

Experimental work on bijels (albeit with larger particles, possibly having a primary at- 
tractive minimum in the pair potential) show these structures to be mechanically arrested on 
time scales of weeks [15] - far beyond the timescales currently accessible in our simulations. 
It is possible that the residual dynamics we observed is the consequence of numerical artefact, 
causing particle ejections when none should be seen; or it may be that the barrier to particle 
detachment from the interfacial layer is roughly two orders of magnitude smaller than was 
expected (a < 0.01 rather than a ~ 1). However, we have argued that the continuing ejec- 
tions and rearrangement seen in our simulations may represent instead a different physical 
effect at intermediate times. They might be caused by a departure from local equilibrium 
during the sequestration and arrest of the colloids, which takes place on a time scale short 
compared to their Brownian time tb, and could lead to a transient population of weakly 
bound particles at the interface which do not all leave until t ^> r#, in a time a regime we 
are yet to reach. 

The work presented here remains consistent with the predicted formation of completely 
arrested structures, including both bijels and frozen droplet phases, even when the colloidal 
interactions are purely repulsive [10]. For this to be sustained, ae must greatly exceed fc#T at 
late times. If the final a value remains of order 0.01 rather than 0(1) as originally expected, 
the minimum particle size for a stable repulsive bijel is increased from a few nanometres to 
ten times this - which is still much smaller than in the experiments done so far [15J. On the 
other hand, we cannot yet exclude the possibility that a is strictly zero, corresponding to 
a crowded monolayer that can continuously shed particles without crossing barriers. This 
would preclude formation of stable bijels in the purely repulsive case. The experimental 
existence of such repulsive bijels thus remains an important open question of principle, 
although from an applications perspective it may be less crucial. After all, strategies to 
engineer bonding interactions between colloids are well developed. Indeed, if repulsive bijels 
are ultimately shown to be unstable, attractive interactions of some form must presumably 
already be present in the experimental system used successfully to make bijels [15J. 
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